library(SpatialEpi)
## Loading required package: sp
library(readxl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
IMR1=read_excel("~/2018 Reported Infant Death.xlsx",skip = 1)
Malesex=read_excel("~/2018 Reported Infant Death.xlsx",sheet = "Gender Male")
head(IMR1)
## # A tibble: 6 × 5
##   Districts   `Place of recorded death` Cases Population (Total Infant B…¹ Age  
##   <chr>       <chr>                     <dbl>                        <dbl> <chr>
## 1 Gaborone    Gaborone                     57                         5986 >28 …
## 2 Gaborone    Gaborone                     23                         2416 1mon…
## 3 Francistown Francistown                  28                         3678 >28 …
## 4 Francistown Francistown                  19                         2495 1mon…
## 5 Lobatse     Lobatse                       8                          496 >28 …
## 6 Lobatse     Lobatse                      10                          621 1mon…
## # … with abbreviated variable name ¹​`Population (Total Infant Born)`
d=aggregate(x=IMR1$Cases,by=list(Districts=IMR1$Districts),FUN=sum)
M=aggregate(x=Malesex$`Gender (M)`,by=list(Districts=Malesex$Districts),
            FUN=sum)
head(d)
##     Districts   x
## 1     Central 303
## 2       Chobe   6
## 3 Francistown  47
## 4    Gaborone  80
## 5      Ghanzi  22
## 6     Jwaneng   6
names(d)=c("id","Y")
names(M)=c("Districts","Male_sex")
pop=IMR1$`Population (Total Infant Born)`
cases=IMR1$Cases
n.strata=3.25
E=expected(pop,cases,n.strata)
d$E=E[match(d$id,unique(IMR1$Districts))]
d=merge(d,M,by.x="id",by.y="Districts")
d$SMR=d$Y/d$E
head(d)
##            id   Y         E Male_sex        SMR
## 1     Central 303  24.67221      158 12.2810226
## 2       Chobe   6  48.55638        3  0.1235677
## 3 Francistown  47  52.64422       25  0.8927855
## 4    Gaborone  80 184.15797       42  0.4344097
## 5      Ghanzi  22  18.08976       12  1.2161576
## 6     Jwaneng   6  26.18266        4  0.2291593
Districts = raster::getData("GADM", country = "Botswana", level = 1)
## Warning in raster::getData("GADM", country = "Botswana", level = 1): getData will be removed in a future version of raster
## . Please use the geodata package instead
library(sp)
rown=sapply(slot(Districts, "polygons"), function(x) slot(x, "ID"))
rownames(d)=rown
identical(rownames(d), getSpPPolygonsIDSlots(Districts))
## Warning: use *apply and slot directly
## [1] TRUE
map= SpatialPolygonsDataFrame(Districts,d,match.ID = TRUE)

head(map@data)
##             id   Y         E Male_sex        SMR
## 1      Central 303  24.67221      158 12.2810226
## 9        Chobe   6  48.55638        3  0.1235677
## 10 Francistown  47  52.64422       25  0.8927855
## 11    Gaborone  80 184.15797       42  0.4344097
## 12      Ghanzi  22  18.08976       12  1.2161576
## 13     Jwaneng   6  26.18266        4  0.2291593
library(leaflet)
l=leaflet(map)%>%
  addTiles()
pal=colorNumeric(palette="YlOrRd",domain=map$SMR)
l %>% addPolygons(color="grey",weight=1,fillColor=~pal(SMR),
                  fillOpacity=0.5) %>% 
  addLegend(pal=pal,values=~SMR,opacity=0.5,title="SMR",
            position="bottomright")
labels=sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/> 
               Male sex proportion: %s <br/> SMR: %s",
               map$id,map$Y,round(map$E,2),map$Male_sex,round(map$SMR,2)) %>%
  lapply(htmltools::HTML)
l %>% addPolygons(color="grey",weight=1,fillColor=~pal(SMR),fillOpacity=0.5,
                  highlightOptions=highlightOptions(weight=4),label=labels,
                  labelOptions=labelOptions(style=list("font-weight"="normal",
                                                       padding="3px 8px"),
                                            textsize="15px",
                                            direction="auto")) %>%
  addLegend(pal=pal,values=~SMR,opacity=0.5,title="SMR",
            position="bottomright")
library(spdep)
## Loading required package: spData
## Loading required package: sf
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(INLA)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## 
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## 
## Loading required package: parallel
## This is INLA_22.12.16 built 2022-12-23 13:24:10 UTC.
##  - See www.r-inla.org/contact-us for how to get help.
nb=poly2nb(map)
head(nb)
## [[1]]
## [1]  2  5  8  9 11 12 13 16
## 
## [[2]]
## [1]  1 12
## 
## [[3]]
## [1] 11
## 
## [[4]]
## [1]  8  9 14
## 
## [[5]]
## [1]  1  7  9 12
## 
## [[6]]
## [1] 15
nb2INLA("map.adj",nb)
g=inla.read.graph("map.adj")
map$re_u=1:nrow(map@data)
map$re_v=1:nrow(map@data)
formula=Y~Male_sex+f(re_u,model="besag",graph=g)+f(re_v,model="iid")
res=inla(formula,family="poisson",data=map@data,E=E,
         control.predictor=list(compute=T))
summary(res)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
##    working.directory = working.directory, ", " silent = silent, inla.mode 
##    = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
##    .parent.frame)") 
## Time used:
##     Pre = 0.854, Running = 0.255, Post = 0.15, Total = 1.26 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.776 0.207     -1.190   -0.775     -0.369 -0.772   0
## Male_sex     0.020 0.004      0.012    0.020      0.029  0.020   0
## 
## Random effects:
##   Name     Model
##     re_u Besags ICAR model
##    re_v IID model
## 
## Model hyperparameters:
##                        mean       sd 0.025quant 0.5quant 0.975quant    mode
## Precision for re_u 18852.32 18808.97    1290.87 13180.04   69010.69 3551.04
## Precision for re_v     3.03     1.26       1.21     2.82       6.10    2.42
## 
## Marginal log-Likelihood:  -93.87 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
head(res$summary.fitted.values)
##                           mean         sd  0.025quant   0.5quant 0.975quant
## fitted.Predictor.01 12.3202871 0.70826119 10.98960406 12.2999794 13.7666008
## fitted.Predictor.02  0.1859291 0.05803335  0.09468564  0.1784237  0.3202189
## fitted.Predictor.03  0.8854982 0.12658889  0.66326213  0.8765663  1.1586771
## fitted.Predictor.04  0.4495950 0.04902685  0.36097573  0.4469962  0.5530276
## fitted.Predictor.05  1.1133586 0.23728791  0.72087100  1.0886917  1.6468377
## fitted.Predictor.06  0.2946768 0.09531422  0.14812262  0.2813314  0.5180600
##                           mode
## fitted.Predictor.01 12.2594222
## fitted.Predictor.02  0.1640752
## fitted.Predictor.03  0.8589968
## fitted.Predictor.04  0.4418447
## fitted.Predictor.05  1.0410090
## fitted.Predictor.06  0.2564490
map$RR=res$summary.fitted.values[,"mean"]
map$LL=res$summary.fitted.values[,"0.025quant"]
map$UL=res$summary.fitted.values[,"0.975quant"]
pal=colorNumeric(palette = "YlOrRd",domain=map$RR)
labels=sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/>
Male sex proportion: %s <br/> SMR: %s <br/> RR: %s (%s, %s)",
               map$id, map$Y, round(map$E, 2), map$Male_sex, round(map$SMR, 2),
               round(map$RR, 2), round(map$LL, 2), round(map$UL, 2)) %>%
  lapply(htmltools::HTML)
leaflet(map) %>% addTiles() %>%
  addPolygons(color = "grey", weight = 1, fillColor = ~pal(RR), fillOpacity = 0.5,
              highlightOptions = highlightOptions(weight = 4), label = labels,
              labelOptions = labelOptions(style = list("font-weight" = "normal",
                                                       padding = "3px 8px"),
                                          textsize = "15px", direction = "auto")) %>%
  addLegend(pal = pal, values = ~RR, opacity = 0.5, title = "RR",
            position = "bottomright")
pal=colorNumeric(palette = "YlOrRd", domain = map$SMR)
leaflet(map) %>% addTiles() %>%
  addPolygons(color = "grey", weight = 1, fillColor = ~pal(SMR), fillOpacity = 0.5,
              highlightOptions = highlightOptions(weight = 4), label = labels,
              labelOptions = labelOptions(style = list("font-weight" = "normal",
                                                       padding = "3px 8px"),
                                          textsize = "15px", direction = "auto")) %>%
  
  addLegend(pal = pal, values = ~SMR, opacity = 0.5, title = "SMR",
            position = "bottomright")